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6 ■ ABSTRACT 

As dust settles in a protoplanetary disk, a vertical shear develops because the 
dust-rich gas in the midplane orbits at a rate closer to true Keplerian than the 
^ \ slower-moving dust-depleted gas above and below. A classical analysis (neglect- 

ing the Coriolis force and differential rotation) predicts that Kelvin-Helmholtz 



instability occurs when the Richardson number of the stratified shear flow is 
below roughly one-quarter. However, earlier numerical studies showed that the 
Coriolis force makes layers more unstable, whereas horizontal shear may stabi- 

(^ \ lize the layers. Simulations with a 3D spectral code were used to investigate 

these opposing influences on the instability in order to resolve whether such lay- 

k>( ; ers can ever reach the dense enough conditions for the onset of gravitational 

H ' instability. I confirm that the Coriolis force, in the absence of radial shear, does 

indeed make dust layers more unstable, however the instability sets in at high 
spatial wavenumber for thicker layers. When radial shear is introduced, the on- 
set of instability depends on the amplitude of perturbations: small amplitude 
perturbations are sheared to high wavenumber where further growth is damped; 
whereas larger amplitude perturbations grow to magnitudes that disrupt the dust 
layer. However, this critical amplitude decreases sharply for thinner, more un- 
stable layers. In 3D simulations of unstable layers, turbulence mixes the dust 
and gas, creating thicker, more stable layers. I find that layers with minimum 
Richardson numbers in the approximate range 0.2 - 0.4 are stable in simulations 
with horizontal shear. 
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INTRODUCTION 



It is a remarkable fact that planets start out as microscopic grains within the proto- 
planetary disks of gas and dust in orbit around newly-formed protostars, somehow growing 
roughly 10^° orde rs of magnitud e in mass in a period no more than 10^ years corresponding 
to disk lifetimes (JLissauerlll993l ). There is no one physical process that can explain growth 
over this enormous range of sizes: the very smallest grains (micron to millimeter sizes) can 
grow via collisional agglomeration in which the sticking mechanism is electrostatic in nature; 
whereas, on the other end of the size spectrum, objects in the kilomet er to tens of kilometers 
regime can grow via gravity-enhanced collisions (JBeckwith et al.ll200Cll ). The least understood 
stage of growth is how millimeter-size particles grow to kilometer-size; grains in this regime 
are too large for sticking via electrostatic forces, yet far too small to have any significant 
self-gravity. Even more problematic is the fact that particles in this intermediate-size regime 
are strongly affected by aerodynamic drag of the surrounding gas: meter-size objects, for 
example, have radial drift speeds on the order of 10'* cm/ s at 1 AU and thus sp iral inward 
onto the protostar on the timescale of a few hundred years (IWeidenschilling||l977l ) . Whatever 
process is responsible for grain growth through this range of sizes must act on timescales 
faster than this inspiral time if any raw materials are to be available to build protoplanets. 



Goldreich fc WardI (Il973l ) and ISafronovl (119691 ) proposed that a very thin, very dense 
sheet of settled dust in the midplane of the protoplanetary disk might be be gravitation- 
ally unstable; the nonlinear evolution would result in the layer clumping-up directly into 
gravitationally-bound kilometer-size planetesimals on a timescale of order the orbital pe- 
riod. According to this scenario, there is a direct jump from small particles to kilometer-size 
planetesimals without growing slowly through the intermediate sizes which have very short 
orbital decay timescales. As attractive as this mechanism is for planetesimal formation, a 
significant obstacle is turbulence which can stir and mix the dust with the gas and prevent 
the dust layer from s ettling into a thin en ough, dense enough sheet for the gravitational 
instability to operate (IWeidenschillingl ll98Cll ) . 



Even in the absence of any mechanism to drive turbulence, the settling of the dust 
particles into a thin sublayer in an initially laminar midplane would create a vertical shear 
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that might be unstable to Kelvin- Helmholtz instability. Because of a relatively weak outward 
radial pressure gradient, pure gas in a protoplanetary disk orbits the protostar at a rate 
slightly sl ower than true Kep l erian: V^„., = Vk(1 —p ), where r] ~ 10~^ for typical conditions 
at 1 AU (JAdachi et al.lll976l : IWeidenschillingI 119771 ). Pure dust, in the absence of any gas. 



would orbit exactly at the Keplerian rate. One can show that in the limit of perfect dust-gas 
coupling, th e orbital velocity of a mixture of gas and dust is a function of the local dust- 
to-gas ratio (JAdachi et al.lll976l ): as dust settles into the midplane, the dust-rich gas in the 
midplane orbits faster than the dust-deple ted gas above and b elow the midplane . Two - fluid 
(gas and dust) nunaerica l simulations by ICuzzi et al.l (Il993l ) , IChampney et al.l (Il995l ) and 
Dobrovolskis et al.l (Il999l ) suggested that turbulent diffusion would prevent settling of grains 
into thin enough sheets for gravitational instability. 

The onset of Kelvin-Helmholtz instability is determined by a competition between the 
stabilizing effects of stratification and the destabilizing effects of vertical shear. If there 
is no rotation and no horizontal shear, then a necessary (but not s ufficient) condition f or 
instability is the Richardson number criterion (jChandrasekhaii Il960l : iDrazin fc Reidlll98ll ): 



instability: Ri{z) 



N' 



(dU/dz) 



< Ricrit ~ 0.25 for some z. 
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where N is the Brunt- Vaisala frequency (the frequency of buoyant oscillations of a stably- 
stratified medium) and dU/dz is the vertical shear. The critical Richardson number of one- 
quarter corresponds to a state in which the kinetic energy of the vertical shear is sufficient 
to lift the heavier gas out of the gravitationally potential well and remix it with the lighter 
overlying fluid. Assuming that dust would set tle into thinn e r and thinner lay e rs dow n to 
the limit set by the c lassic Richardson criterion. ISekiyal (ll998l ). ISekiya fc Ishitsul (|2000| ). and 
Youdin &: Shij (120021 ) determined critical vertical quasi-equilibrium dust profi les and investi- 
gated the conditions necessary for such layers to be gravitationally unstable. iGaraud &: Lin 
( 120041 ) pursued two-fiuid linear calculations (without rotation and radial shear) of dust sed- 
imentation into sheets and found that the Kelvin-Helmholtz instability was excited before 
gravitation al instability unless the g lobal dust-to-gas ratio was greatly enhanced over solar 
abundance. lYoudin fc Chiang] (12004 ) suggested that such enhancements might be attainable 
because the inward drift speed decreases as particles migrate inward, resulting in "particle 
pile-ups." 

It is not at all clear whether the classic Richardson number criterion is appropriate for 
the Kelvin-Helmholtz instability in protoplanetary disks in which both rotation and radial 
shear could significantly affect t he onset of instabi l ity as well as alter the nonlinear evolution 
of any instability that develops. Ilshitsu fc Sekiyal (120031 ) pursued a purely linear analysis to 
study the effect of radial shear on the time evolution of unstable Kelvin-Helmholtz modes 
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and showed that the modes were s heared to higher spa t ial wa venumber and could eventually 



be stabilized. On the other hand, iGomez &: Ostrikerl (120051 ) investigated the effects of the 



Coriolis force on the onset of instability, and found that settled layers were unstable at much 
higher Richardson numbers (corresponding to thicker dust layers) than in the classic case. 
However, these simulations were two-dimensional and did not include the radial shear. 

The motivation for this work is to investigate the apparently opposing influences of 
rotation (destabilizing) and radial shear (stabilizing) on the Kelvin-Helmholtz instability of 
settled dust layers and resolve whether or not such layers can ever reach the thin enough, 
dense enough conditions for the onset of gravitational instability. Apart from the impact of 
Kelvin-Helmholtz instability on the formation of planetesimals, it is importan t to determine 



the yertical distribution of dust in order to co rrectly interpret observations (JBrittain et al. 



20051 : iRettig et al.ll2006l : iDuUemond et al.l 120071 ). The approach we take here is that the dust 



is perfectly coupled to the gas; in terms of timescales, the friction or stopping time (see 
(12-81) below) is much shorter than the evolution of Kelvin-Helmholtz instabilities if they are 
present. We will not simulate the formation of the dust layer from a well-mixed state, but 
assume that the layers have already formed with a given profile of the dust-to-gas ratio. If 
the layer turns out to be unstable on a fast timescale, it imphes that such a layer would never 
have formed in the first place. In §2, we present equations for the hydrodynamic evolution of 
settled dust layers in the limit of perfect dust-gas coupling. In §3, we revisit the cases of (i) 
no rotation, no shear, and (ii) rotation, no shear. In §4, we present new three-dimensional 
simulations of the evolution of settled dust layers with both rotation and shear. Finally, in 
§5, we discuss the impact of these results on the planetesimal formation via gravitational 
instability. 



2. SINGLE-FLUID EQUATIONS FOR PERFECTLY COUPLED GAS & 
DUST IN 3D CARTESIAN SHEARING BOX 

2.1. Equilibrium for a Gas Disk 

Consider the time- in dependent, axisymmetric azimuthal flow V^ of gas around a proto- 
star of mass M^: 

V^ 1 

^f = V<l> + —VPg, (2-1) 

r Pg 

where (r, 0, z) are protostar-centered cylindrical coordinates with corresponding unit vectors 
(f , 0, z), $ = —GMi,j\/r'^ + z^ is the gravitational potential, G is the gravitational constant, 
and pg and pg are the equilibrium gas density and pressure. Protoplanetary disks are ther- 
mally cool in the sense that the gas sound speed Cg is much slower than the Keplerian orbital 
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velocity Vxir) = rrixir) = ^jGM^Jr. Hydrostatic balance implies that the time it takes 
sound waves to traverse the thi ckness of the disk is of order the orbital period; thus, cool 



disks are geometrically thin (see lFrank et al.lll985l ): 



c, ~ ^Kiig, (2-2a) 

6 = c,IVk ~ iijr<\, (2-2b) 

where Hg is the vertical pressure scale height. In cool disks, the radial component of the 
protostellar gravity nearly balances the centrifugal force, but because of the relatively weak 
outward radial pressure force, the gas orb its at slightly slower than the Keplerian velocity 



flAdachi et all Il976l : IWeidenschillind 119771 ): 



fi(r,z) = fi^(r)[l-r7(r,z)], (2-3a) 



-{dpg/dr)/pg ^ ^ fz 



2 



where the fractional deviation from Keplerian is r/ ~ 5^ ^ 1. Depending on the disk model, 
7] typically takes values between 10^'^ — 10~^ at r = 1 AU, and the Mach number for the 
maximum deviation from Keplerian is of order vjVk/cs ~ 0.1 

In contrast to a pure gas disk which orbits at a slightly sub-Keplerian speed, a disk 
of pure ( "pressureless" ) dust orbits at the full Keplerian speed. As dust sediments into a 
thin sheet and creates a vertical shear, we expect that Kelvin-Helmholtz instabilities might 
potentially develop when the local dust-to-gas ratio /i = Pd/pg in the midplane approaches 
of order unity: 

/io = M^ = 0)~l when X = Hd/Hg^ Ed/Eg, (2-4) 

where pd is the local mass density of dust, Hd is the scale height of the dust sub-layer, T^d 
is the surface mass density of dust and S^ is the surface mass density of gas. When Kelvin- 
Helmholtz instabilities do develop, the horizontal and vertical scales of the fastest growing 
unstable modes are typically of order the thickness of the shear layer: (Ar, r A0, Az) ~ 
Hd ~ ^Hg. We can divide the velocity into two components: the velocity across the region 
of interest due to the Keplerian shear v ~ Hd^K ~ XHgQ ~ Xcg, and the differential velocity 
between the settled dust sub-layer and the dust-depleted gas above and below the midplane 
V ~ tjVk ~ Scg. The Mach number of the flow is thus: e = max(5. A). 

Motivated by this dimensional analysis, we simulate the dynamics only within a small 
patch of the disk (r — tq, ro(0 — 0o), z) -^ (x, y, z) ^^ Hd ^ XHg that co-rotates with the gas 
at some fiducial radius rg with angular speed Qp = fixoi^ ~ Vo): where Qko = ^k(^o) and 
rjQ = ri{rQ, z = 0). The tidal term (i.e., the remainder after the near cancellation of the inward 
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radial protostellar gravity and the outward centrifugal force) and the equihbrium pressure 
gradient are given by: 

-^ + nlr = nl,r,[3x/ro-2r]o + 0{5^\\5^\,5^)], (2-5a) 

-y^ = nl,r42r^o + 0{6'X',6'X)], (2-5b) 

--f^ = ^loro[z/ro + 0{6'X')]. (2-5c) 

In this work, we assume the background gas temperature is spatially constant, T = Tq, so 
that the gas pressure gradient can be written in terms of the gas density gradient: (Vpg)/pg = 
TZTqV In Pg, where TZ = Cp — Cy is the gas constant. The equilibrium gas density is thus: 

Pg{x,z) = pQexp[-x/Ag- z^/2Hg) , (2-6a) 

H' ^ 7^To/fi^o, (2-6b) 



Ag = nTo/2nj,Q7]oro = {cso/2r]oVKo)Hg. (2-6c) 



-PgV ■ V, (2-7b) 

-PdV ■ V, (2-7c) 



2.2. Dynamic Equations for Gas &; Dust 

The Euler equations for perfectly coupled gas and dust in the 3D Cartesian shearing 
box are: 

-— = -2nKoz X V + {snl^QX - 2n'l^Qror]o) x - Vi'^kqZz - ■ rVpg, (2-7a) 

at [Pg + pd) 

dpg 

dt 

dpd 

dt 

r/s 

_^ = 0, Sg^CyHpgp-r) + Sgo (2-7d) 

where v is the velocity of a parcel of gas and dust, Pg, Pg, and Sg are the pressure, density, 
and entropy of the gas; pg is the density of the dust; 7 = Cp/Cy is the ratio of specific heats 
at constant pressure and constant volume; the advective or Lagrangian derivative is defined 
d/dt = {d/dt + V ■ V). The key difference between the dynamics of the gas and that of the 
dust is that we treat the dust as a cold, pressureless fluid. 

The dust continuity equation can be recast in terms of the local dust-to-gas ratio: 
dp/dt = 0, that is, the local dust-to-gas ratio is an advectively conserved quantity, meaning 
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that a parcel of fluid maintains its dust content. This is simply a consequence of the perfect 
coupling assumption which does not allow the dust component to slip apart from the gas 
component. This approximation is valid if the stopping time ts (the e-folding time for a 
particle's velocity to match that of the surrounding medium because of frictional coupling) 
is much shorter than the timescales of interest for Kelvin-Helmholtz instability. For small 
particles for which the gas mean-free-path is larger than the size of the particles, Epstein 
drag sets this timescale (ICuzzi et al.lll993l : iGaraud et al.l 120041 ): 



t 



S,epstein/ l^orb 



/to 



27r/fi 



10" 



K 



1 cm 



r \3/2 
1 AuJ 



(2- 



where ps and a are the solid density and radius of a dust grain, and pg and Cg are the gas 
density and sound speed. 

Because we expect the Kelvin-Helmholtz instability to set in at low Mach number e ~ 
max((5. A) ^ 1, we invoke the anelastic approximation for the gas flow. We decompose the 
gas pressure and density into their equilibrium components (denoted with overbars) and 
fluctuating components (denoted with tildes): Pg = pg + pg, Pg = Pg + Pg- At low Mach 
number, the fluctuating components should scale as: pg/pg ~ Pg/Pg ~ ^^ << 1- The gas 
pressure gradient and gas continuity equations can then be expanded: 







— ypg + —'^Pg-——^Pg 

Pg Pa P9P9 

[V ■ (p,v)] [1 + 0{^ 



[l + 0{t 



(2-9a) 
(2-9b) 



The anelastic approximation has been used extensively in the study of deep, s ubsonic convec- 

tion in planetary atmospheres (Qgura fc Phillips 19621: lGoughlll969) and stars (Gilman &: Glatzmaier 



J2000l . 



1981 



2005 



Glatzmaier &: Gilmanlll981a| ]bl) . iBarranco et al.l (120001 ) and iBarranco fc Marcus 



20061 ) previously used the anelastic approximation to study 3D vortices in protoplane- 



tary disks. One of the consequences of this approximation is that the total density is replaced 
by the time-independent mean density in the mass continuity equation, which has the effect 
of filtering high-frequency acoustic waves and shocks, but allowing slower wave phenomena 
such as internal gravity waves. 

The dynamic equations for coupled gas and dust with the constant temperature back- 
ground and anelastic approximation become: 



T rji 

-— = -2nKoz X V + 31l|-oXx + — (211^oror7ox + ^ko^z) - Vhg 
dt Iq 



p 



p+l 



1 + — I {2QJ^Qror]oit + ^koZz) - Vhg 



(2-lOa) 



= V-p<,v, 

/i = Pd/pg, 



^ = 

dt 
df 



Pg 



1 + — I V ■ (2n^oro?7ox + n^Ko^z) 



PgTZTo + PglZf, hg = Pg/Pg 



(2-lOb) 
(2-lOc) 

(2-lOd) 
(2-lOe) 



When the background temperature is spatially constant, the gas enthalpy turns out to be 
be a more useful quantity than the gas pressure: hn = P n /Pg- These equations are almost 
identical to the anelastic equations in iBarranco fc Marcus! (J2005|, l2006l ) , with the addition of 
a nonlinear forcing term for the inertia of the dust. 

The above set of dynamic equations allow the following steady-state equilibrium (de- 
noted with the dagger symbol): 





/.t(z) 

vl{x,z) 

dhliz) 
dz 



,t 



f\ 



arbitrary, 



■|r2ii-ox + 






^KoroVo, 



-fi^{z)nj^oZ. 



(2-lla) 
(2-llb) 

(2-llc) 
(2-lld) 



In the limit where we take Ag-^oo and neglect the radial variation of the background 
gas density and radial component of the gas buoyancy, one can derive the following global 
energy balance equations: 



KE 
PE 



- {KE + PE) 



f l{l + ^l)pg^r■^rdV, 
Jv 

yufi^o (2^0^70^ + 1^;^ - fa;^) + Cppgf 

[(1 + p)pgV^Vy] \x=+L^/2 dydz 

2Q^j^Qror]oLx / [pgpv^] \x=+l^/2 dydz, 
Js 



dV, 



-qP-koLx 



(2-12a) 
(2-12b) 

(2-12c) 
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2.3. Brief description of numerical method 

Her e we briefly describe the n umerical method; a more detailed presentation can be 
found in lBarranco fc Marcus! (120061 ) . We solve the dynamic equations fl2-10p with a spectral 
method; that is, eac h variable is represented as a flnite surn of basis functions rnultiplied b y 
spectral coeflicients (JGottlieb k. Orszagill977l : lMarcuslll986l : ICanuto et al.lll988l : lBoydlll989l ). 
The choice of basis functions for each direction is guided by the corresponding boundary 
conditions. The equations are autonomous in the azimuthal coordinate y, so it is rea- 
sonable to assume periodic boundary conditions in this direction. However, the equations 
explicitly depend on the radial coordinate x because of the linear background shear. We 
adopt "shearing box" boundary conditions: q{x + Lx,y — {3/2)flQLxt,z,t) = q{x,y,z,t), 
where q represents any of v, /z., T, etc. In practice, we rewrite the equations fl2-10p in 
terms of quasi-Lagrangian or shearing coordinates that advect with th e background shear 
JColdreich fc Lvnden-Bell Il965l : iMarcus fc Pre3 Il977l : iRogallol Il98lh : f = t, x' = x, 
y' = y + {3/2)flQxt, and z' = z. In these new coordinates, the radial boundary condi- 
tions become: q{x' + L^, y', z', t') = q{x', y', z', t'). That is, shearing box boundary conditions 
are equivalent to periodic boundary conditions in the shearing coordinates. Physically, this 
means that the periodic images at different radii are not flxed, but advect with the back- 
ground shear. 



In the shearing coordinates, the equations are autonomous in both x' and y' (although 
they now explicitly depend on t'), making a Fourier basis the natural choice for the spectral 
expansions in the horizontal directions: 



q{x',y',z',t') = J2Ut')e'''^'''e"''y''U^'), 



(2-13) 



where q is any variable of interest, {^k(^')} is the set of spectral coeflicients, and k = 
{k'^, k' n} is the set of wavenumbers. We have implemented the simulations with two differ- 
ent sets of basis functions for the spectral expansions in the vertical direction, corresponding 
to two different sets of boundary conditions: 



(i) For the truncated domain —L^ < z < L^ 

Tn{z/Lz) = cos(n^), where ^ = cos 



-ii 



z/L^ 



we use Chebyshev polynomials: (f)n{z) = 
We apply the condition that the vertical 



velocity vanish at the top and bottom boundaries: Vz{x,y,z = ±Lz,t) = 0. 

(ii) For the inflnite domain — oo < 2; < oo, we use rational Chebyshev functions: (f)n{z) 



cos(n^) (for v^, Vy, and all the thermodynamic variables) or 



sin(n^) (for Vz), where 



^ = cot~^{z/Lz). In this context, Lz is no longer the physical size of the box, but is a 
mapping parameter; exactly one half of the grid points are within \z\ < Lz, whereas the 
other half are widely spaced in the region Lz < \z\ < cxd. No explicit boundary conditions 
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on the vertical velocity are necessary when we solve the equations on the infinite domain 
because the (pn{z) = sm{n^) basis functions individually decay to zero at large z. 

The equations are integrated forward in time with a fractional step method: the nonlin- 
ear advection terms are integrated with an explicit second-order Adams-Bashforth method, 
and the pressure step is computed with a semi-implicit second-order Crank-Nicholson method. 
The time integration scheme is overall globally second-order accurate. Unlike finite-difference 
methods, spectral methods have no inherent grid dissipation; energy cascades to smaller and 
smaller size scales via the nonlinear interactions, where it can "pile-up" and potentially 
degrade the convergence of the spectral expansions. We employ a V^^ hyperviscosity or 
low-pass filter every timestep to damp the energy at the highest wavenumbers. 

Different horizontal Fourier modes interact only through the nonlinear advective terms; 
once these terms are computed, the horizontal Fourier modes can be decoupled. This mo- 
tivated us to parallelize the code: each processor computes on a different block of data in 
horizontal Fourier wavenumber space. Parallelization is implemented with Message Passing 
Interface (MPI), typically using between 64 and 512 processors. Wall-clock time scales in- 
versely with number of processor s , indi cating near-optimal parallelism; timing analyses are 
presented in iBarranco &: Marcus! (120061 ) . 



3. REVISITING THE CASE OF NO DIFFERENTIAL ROTATION 



The linear stability of settled dust layers to Ke lvin-He l mholtz instability has previ 



2000 



ously been investigated by a number of researchers (ISekival Il998l : ISekiya fc Ishitsu 
Youdin fc Shull2002l : iGaraud fc LinI |2004J : iGomez fc Ostrikeii l2005l ) : however, almost all of 



these analyses neglected the role of the differential rotation in the protoplanetary disk. If 
one were to directly linearize the equations (12-lOp . one would find the resulting set to de- 
pend linearly on the radial coordinate x because of the shear, making it difficult to apply 
periodic boundary conditions in the radial direction. Alt ernatively, one could employ a 
set of coordinates that advect with the background shear (JGoldreich &: Lynden-Belll Il965 



Marcus fc Press! 119771 : !Ryu fc Goodman 



19921). but t hen th e resulting linearized equations 



would explicitly depend on time. !lshitsu fc Sekiya! (12003! ) employed this approach using 
an initial-value code to address the effect of horizontal shear on unstable eigenmodes in 
the small-amplitude, linear regime. They found that modes grew for a limited period of 
time, but were eventually stabilized as they were sheared out to high spatial wavenumber. 
No matter the approach, the linear stability analysis is difficult to treat analytically, mo- 
tivating theorists to tackle the problem without the background shear with the hope that 
the results are qualitatively, if not quantitatively, useful in determining whether or not the 
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Kelvin-Helmholtz instability is a barrier to further settling. 

In this section, we briefly revisit the case of no radial shear, flrst without, and then 
with the Coriolis force. We will neglect the x (radial) dependence of the background gas 
density (setting A^ — » oo), which eliminates the radial component of gas buoyancy. We 
assume eigenmodes of the form g'(t, x, y, z) = q\z) exp(— icut + ik^x + ikyy). The linearized 
equations corresponding to (12-lOp . without radial shear and radial gas buoyancy, are: 



-lUV'y 

-iujv'^ 

r 

To 

—iujfi' 





-vl{z)ikyv'^ + 2QKoVy - 
-2Qkov'^ - vl{z)ikyVy - 
-vl{z)ikyv'^ 



%Kx 



[l + ^^Kz)] 



H 



h' 



{d/dz) 



T' 

-vliz)iky— + 

-vl{z)ikyfi' - 



TZ dlnpg 
Cp dz 
dfi'' 
dz 



dz 
h' + 



V, - 



tfvi 



r 



2Ql^Qr]oro 

[l + /it(^)p 

h\ 



/^ 



02 7 



L[i + /^n^)] 



02 7 



To L[l + /it(z)]J 



d dlnpg 
dz dz 



(3-la) 

(3-lb) 

^'(3-lc) 
(3-ld) 
(3-le) 
(3-lf) 



One can eliminate T' and fi' from the equation for v'^, yielding: 



[-iu + vl{z)iky]\'^ 



-iuj + vl{z)iky 



dh 



+ 



Oj^qZ 



[1 + ^n^)] 
This form allows us to easily identify the Brunt- Vaisala frequency 

N^ -- 



dz ' [l + fi^z)] \Cp dz dz ' '' ^ ' 



02 7 



IZ dlnpg dp^ 
[1 + p^{z)] \C^ dz ^ H 



(3-3) 



Before proceeding with the linear analysis, we must specify the form of the vertical 
distribution of dust. We choose a Gaussian proflle for the local dust-to-gas ratio: 



/xt(z) 



plexp{-zy2Hl), l^l^^^, H-, 



Hf - H-\ (3-4) 



where we have deflned the midplane dust-to-gas ratio /Iq, the initial Gaussian scale height 
for the dust density if^, and the initial Gaussian scale height for the dust-to-gas ratio H^. 
The gradient Richardson number (ll-lj) for a Gaussian distribution of dust is: 



Riiz) 



VoVi 



KO 



CsO 






p^{z) 



n 1 



Cpp'^{z) \H, 



H,, 



(3-5) 
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Fig. 1. — Gradient Richardson number (13-51) for Gaussian profiles of tlie local dust-to-gas ra- 
tio /i. For these three curves, S^/E^ = 0.01, tiqVko/cso = 0.1, and H^j Hg = 0.04, 0.02, 0.005, 
corresponding to peak dust-to-gas-ratios of /Iq = 0.25, 0.50, 2.0. For /Iq = 1/2, the Richard- 
son number is at a minimum right at the midplane {z = 0) and is nearly constant throughout 
much of the dust layer before rising sharply at the edge of the dust distribution. For /xj > 1/2, 
the Richardson number is still relatively constant throughout the core of the layer, although 
the minimum has shifted off the midplane. 
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In Figure [T|, we have graphed the gradient Richardson number with Hd/^g = 0.01 and 
tjoVko/cso = 0.1, for a few different dust scale heights. For yuj = 1/2 {Hd/Hg = 0.02), the 
Richardson number is at a minimum right at the midplane {z = 0) and is nearly constant 
throughout much of the dust layer before rising sharply at the edge of the dust distribution. 
For nl > 1/2 {Hd/Hg < 0.02), the Richardson number is still relatively constant throughout 
the core of the layer, although the minimum has shifted off the midplane. For very thin dust 
layers, we can ignore the gas buoyancy near the midplane (i.e., the TZ/Cp term on the far 
right of equation fl3-5p ): the minimum gradient Richardson number is thus: 



RL 



VoVi 



KO 



CsO 



En 



^^ (27/4) at z = ±H^^2\n{2^^l) for /xj > 1/2 ^^_^.^ 
[I + 4)74 at 2 = for 4 < 1/2. 



It is interesting to note that for dust-rich layers (/ig > 1/2), the minimum Richardson number 
is independent of dust-to-gas ratio. 



We numerically solve the eigenproblem (13- ip for the complex frequencies uj with a 
Cheby shev spectral method in which the top and bottom bo undaries are mapped to in- 
finity (JBarranco fc Marcus! l2006l : iBoydl Il989l : ICain et al.l Il984l ) . Because we are neglecting 
differential rotation, we invoke Squire's theorem which states that two-dimensional eigen 



modes are more uns table than three-dimensional ones (jSquird Il933l : I Chandr asekharl Il961 
Drazin fc Reidlll98ll ). and so set kx = 0. We also restrict our analysis to modes for which 



the dust-to-gas ratio perturbation /i' is an odd-function of the vertical coordinate z. 



3.1. Case of no Coriolis force and no horizontal shear 



Figure [2] shows the growth rates (imaginary part of the complex eigenvalue u) for the 
case where the Coriolis force is turned-off. Each of the twenty plots corresponds to different 
values of the global dust-to-gas ratio, S^/S^, and the strength of the radial gas pressure 
gradient, tiqVko/cso, which sets the maximum differential velocity between pure dust and 
pure gas . Global dust-to-gas ratio varies across rows with values S^/E^ = 0.08, 0.04, 0.02, 
0.01, 0.005. Strength of radial pressure gradient varies down columns with values ?7o^o/cso 
= 0.2, 0.1, 0.05, 0.025. The horizontal axis of each individual plot is the nondimensionalized 
wavenumber (8/7r) iJ^/cy, and the vertical axis is the ratio of the dust scale height to the gas 
scale height. The solid contours, from outer to inner, correspond to growth rates of 0.1, 
0.2, 0.3, 0.4, 0.5, 0.6 in units of f2]^o- "^^^ outermost (dotted) contour corresponds to an 
extrapolation to zero growth rate. 

The peak of the zero-growth contour reveals the thinnest layer to remain stable to 
Kelvin-Helmholtz instability as well as the wavelength of the eigenmode at the onset of 
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Fig. 2. — Contour plots of growth rates of Kelvin-Helmholtz instability for the case of no 
Coriolis force and no horizontal shear. Global dust-to-gas ratio varies across rows with 
values Sd/Sg = 0.08, 0.04, 0.02, 0.01, 0.005. Strength of radial pressure gradient varies 
down columns with values r/oVft-o/cso = 0.2, 0.1, 0.05, 0.025. The horizontal axis of each plot 
is the nondimensionalized wavenumber {S/'K)H(iky, and the vertical axis is the ratio of the 
dust scale height to the gas scale height. The solid contours, from outer to inner, correspond 
to growth rates of 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 in units of I^^^q. The outermost, dotted contour 
corresponds to an extrapolation to zero growth rate. 
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Fig. 3. — Contours in the (Srf/Sg,?7oVft'o/cso) plane for the minimum dust layer thickness 
H(i/ Hg (solid black lines) and minimum Richardson number Rimin (dotted lines) at the onset 
of instability. Over the range of parameter space explored, the minimum Richardson number 
is in the range 0.18 to 0.25, as expected for "classic" Kelvin-Helmholtz instability with no 
Coriolis force and no horizontal shear. Note how the contours for H^/Hg are nearly vertical, 
indicating a weak dependence on the glo bal dust-to-gas ratio for the onset of instability. 
This figure is consistent with Figure 11 in lGaraud fc LinI (120041 ). 
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instability. The wavelength at onset is typically between 8 and 16 times the dust scale 
height over the range of parameter space explored. In Figure [3|, we plot contours in the 
(Srf/Sg, r/o^o/cso) plane for the minimum dust scale height H^/ Hg (solid black lines) to 
remain stable. We also plot contours (dotted lines) for the minimum gradient Richardson 
number (l3-6p corresponding to the minimum dust thickness at the onset of instability. Over 
the range of parameter space explored, the minimum Richardson number is in the range 0.18 
to 0.25, as expected for "classic" Kelvin-Helmholt z instability with no C oriolis force and no 



horizontal shear. Similar results were obtained by iGaraud fc LinI (120041 ) (Figure 11 in their 
work) . 

We simulate the nonlinear evolution of the instability in order to investigate the nature 
of the subsequent mixing of gas and dust. An example of a two-dimensional simulation in 
the y — z plane with S^/Sg = 0.01 and r^oV/co/cso = 0.1 is presented in Figure HI The 
initial dust scale height is Hd/Hg = 0.01, corresponding to a peak local dust-to-gas ratio 
in the midplane of /Iq = 1 and a minimum Richardson number of Rimm = 0.0675. The 
first column illustrates the evolution of the local dust-to-gas ratio /i (deep red = 1, deep 
blue = 0); the second column shows the evolution of the radial component of vorticity 
ujx = dvz/dy — dvy/dz (red = vorticity that points into the page, blue = vorticity that 
points out of the page). The times corresponding to each frame, in units of the orbital 
period, are: 3.8, 4.5, 5.1, 5.7, 15.3. The dust layer develops waves, which grow and break 
into pairs of anti-aligned vortices. Characteristic of two-dimensional turbulence, like-signed 
vortices merge to form larger vortices as energy cascades to larger spatial scales. These 
vortices chaotically interact, leading to thorough mixing of the dust with the gas. 



3.2. Case with Coriolis force, but no horizontal shear 

Figure [5] shows the growth rates for the case with the Coriolis force, but no horizontal 
shear. As before, each of the twenty plots corresponds to different values of the global dust- 
to-gas ratio, Srf/Sg, and the strength of the radial gas pressure gradient, tjoVkq/csq, which 
sets the maximum differential velocity between pure dust and pure gas. Global dust-to-gas 
ratio varies across rows with values S^^/Sg = 0.08, 0.04, 0.02, 0.01, 0.005. Strength of radial 
pressure gradient varies down columns with values r/oVft-o/cso = 0.2, 0.1, 0.05, 0.025. The 
horizontal axis of each individual plot is the nondimensionalized wavenumber {8/iT)Hdky, 
and the vertical axis is the ratio of the dust scale height to the gas scale height. The solid 
contours, from outer to inner, correspond to growth rates of 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 in 
units of ^x\,- Note that the vertical scale of each plot is a factor of 2.2 larger than the 
corresponding ones in Figure [2], and the horizontal scale is a factor of 2 larger. 
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Fig. 4. — Two-dimensional nonlinear evolution of Kelvin-Helmholtz instability with no Cori- 
olis force and no horizontal shear. In this simulation, S^/Sc, = 0.01 and fjoVKo/cso = 0.1. 
The initial dust scale height is Ha/Hg = 0.01, corresponding to a peak local dust-to-gas ratio 
in the midplane of /Xg = 1 and a minimum Richardson number of Rimm = 0.0675. The first 
column illustrates the evolution of the local dust-to-gas ratio /i (deep red = 1, deep blue = 
0); the second column shows the evolution of the radial component of vorticity uj^ (red = 
vorticity that points into the page, blue = vorticity that points out of the page). The times 
corresponding to each frame, in units of the orbital period, are: 3.8, 4.5, 5.1, 5.7, 15.3. 
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Fig. 5. — Contour plots of growth rates of Kelvin-Helmholtz instability for the case with 
the Coriolis force, but no horizontal shear. Global dust-to-gas ratio varies across rows with 
values Eci/Eg = 0.08, 0.04, 0.02, 0.01, 0.005. Strength of radial pressure gradient varies down 
columns with values tiqVko/cso = 0.2, 0.1, 0.05, 0.025. The horizontal axis of each plot is 
the nondimensionalized wavenumber {8/7T)Hdky, and the vertical axis is the ratio of the dust 
scale height to the gas scale height. The solid contours, from outer to inner, correspond 
to growth rates of 0.1, 0.2, 0.3, 0.4, 0.5, 0.6 in units of f^^^g. Note that the horizontal and 
vertical axes of each plot have roughly twice the range as the corresponding plots in Figure [2] 
for the case of no Coriolis force. 
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Gomez &: Ostrikerl (120051 ) were the first to note that settled dust layers are more unstable 
when the Coriolis force is included. Instability occurs for thicker layers and for a much larger 
range of wavenumbers. In the case where there is no Coriolis force, the range of unstable 
wavenumbers for a given ratio of H^/Hg is relatively narrow. In contrast, when the Coriolis 
force is included, we find that instability occurs for very large wavenumbers, with no apparent 
upper limit. However, these large wavenumber (small wavelength) eigenmodes are the ones 
which will be most affected by the inclusion of horizontal shear. 

Figure [6] shows the two-dimensional nonlinear evolution of Kelvin- Helmholtz instability 
with Coriolis force but still no horizontal shear. This simulation is exactly the same as 
the one in Figure HI except that the Coriolis force is included. As before S^/Sg = 0.01, 
tjoVko/cso = 0.1, Hd/Hg = 0.01, /ig = 1 and Rimin = 0.0675. The first column illustrates 
the evolution of the local dust-to-gas ratio /i (deep red = 1, deep blue = 0); the second 
column shows the evolution of the radial component of vorticity Ux (red = vorticity that 
points into the page, blue = vorticity that points out of the page). The times corresponding 
to each frame, in units of the orbital period, are: 1.9, 2.2, 2.5, 2.9, 7.6. Waves appear on the 
dust layer as with the case with no Coriolis force, but no large-scale vortices develop. The 
vorticity has more power at the smallest spatial scales. The nonlinear mixing, however, is 
still very efficient, and the dust is completely re- mixed with the gas throughout the entire 
computational domain. 

We also simulate a thick dust layer that would have been unambiguously stable if there 
was no Coriolis force. Figure [7] shows the two-dimensional nonlinear evolution for the case: 
Srf/Sg = 0.01, tjoVko/cso = 0.1, Hd/Hg = 0.04, /Xq = 0.25 and Rimin = 1.25. The times 
corresponding to each frame, in units of the orbital period, are: 3.2, 3.8, 4.5, 5.1, 7.6. The 
instability sets in at a very small wavelength, fully consistent with the value determined 
in the linear stability analysis in Figure O In the next section, we explore the effect of 
horizontal shear on the evolution of the instability of such thick layers. 



4. 3D SIMULATIONS WITH RADIAL SHEAR 

We present a series of fully three-dimensional simulations of settled dust layers with 
the Coriolis force and differential rotation; Table [1] contains a listing of parameters for these 
simulations. Dust layers in equilibrium were initialized according to (12-111) and (l3-4p . and 
then perturbations were added to the dust-to gas ratio /x: 

/i(x, y, z) =ii\z){l + A{x, y) [cos{nz/2H^) + sm{TTz/2H^)]} . (4-1) 
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Fig. 6. — Two-dimensional nonlinear evolution of Kelvin- Helmholtz instability with Coriolis 
force but still no horizontal shear. In this simulation, S^^/S^ = 0.01 and fjoVKo/cgo = 0.1. 
The initial dust scale height is Ha/Hg = 0.01, corresponding to a peak local dust-to-gas ratio 
in the midplane of /ij = 1 and a minimum Richardson number of Rimm = 0.0675. The first 
column illustrates the evolution of the local dust-to-gas ratio /i (deep red = 1, deep blue = 
0); the second column shows the evolution of the radial component of vorticity uj^ (red = 
vorticity that points into the page, blue = vorticity that points out of the page). The times 
corresponding to each frame, in units of the orbital period, are: 1.9, 2.2, 2.5, 2.9, 7.6. 
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Fig. 7. — Two-dimensional nonlinear evolution of Kelvin- Helmholtz instability with Coriolis 
force but still no horizontal shear. In this simulation, S^^/S^ = 0.01 and fjoVKo/cgQ = 0.1. 
The initial dust scale height is Ha/Hg = 0.04, corresponding to a peak local dust-to-gas ratio 
in the midplane of /ij, = 0.25 and a minimum Richardson number of Rimm = 1-25. The first 
column illustrates the evolution of the local dust-to-gas ratio /i (deep red = 0.25, deep blue 
= 0); the second column shows the evolution of the radial component of vorticity uj^ (red = 
vorticity that points into the page, blue = vorticity that points out of the page). The times 
corresponding to each frame, in units of the orbital period, are: 3.2, 3.8, 4.5, 5.1, 7.6. 
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The amplitude function A{x, y) is constructed in wavenumber space so that each Fourier 
mode has random phase and an amphtude inversely proportional to horizontal wavenumber: 
A{k±) oc fcj^. The perturbations were also forced to be antisymmetric about the x axis 
so that the initial kinetic and potential energies fl2-12p were unchanged. In Table [T], the 
amplitude Arms is the root-mean-square (rms) of these dust-to-gas ratio perturbations in 
the midplane z = 0. 



4.1. The dependence on initial amplitude of perturbations 

Figure [H] shows the results of Run 25, with T^^/T^g = 0.01, tiqVko/cso = 0.1, H^/Hg = 
0.01, yuj = 1, and initial Rimin = 0.0675. These are the same parameters as those used in the 
2D runs with no horizontal shear shown in Figure H] (no Coriolis force) and Figure [6] (with 
Coriolis force). The perturbations in Run 25 had initial amplitude Arms = 0.01 (see also Run 
01 at a lower resolution). Run 26 had exactly the same parameters as Run 25, except the 
amplitude of perturbations was reduce by a factor of 10 (see also Runs 02 and 03). This layer 
is unstable according to the classical Richardson criterion, as demonstrated in the previous 
section for the cases without horizontal shear. However, with the addition of horizontal shear, 
the stability of the dust layer depends also on the amplitude of the initial perturbations: if the 
magnitude of perturbations is below some threshold, the layer remains stable; whereas if the 
amplitud e exceeds some c r itical amount, the layer suffers Kelvin- Helmholtz instability. As 



found by llshitsu fc Sekiyal (120031 ). unstable eigenmodes have a finite period of growth before 
the shear stretches them out to high wavenumber and damps further growth. The nonlinear 
evolution of a dust layer depends on whether the unstable eigenmodes were able to reach 
a sufficient amplitude to trigger nonlinear interactions, resulting eventually in turbulence 
and mixing of the dust with the gas. We have found through experimentation that the 
exact critical amplitude depends on such factors as the kind and shape of perturbations 
(e.g., perturbations to dust-to-gas ratio, or temperature, or velocity field), the resolution 
(number of spectral modes), and the kind and magnitude of small-scale dissipation (i.e., 
hyperviscosity) in the code. 

In Runs 06 and 07, we make the dust layer 50% thicker than the layer in Run 25: 
Sd/Sg = 0.01, VoVko/cso = 0.1, Hd/Hg = 0.015, /ij = 0.667, and initial Rimin = 0.152. The 
amplitudes of initial perturbations were Arms = 0.1 and Arms = 0.04, respectively. A layer 
of this thickness would be unstable according to the Richardson criterion. We again find 
that the stability depends on the amplitude of perturbations. Because this layer is closer to 
stability than the one in Run 25, the eigenmodes would have a slower rate of growth; we 
would expect that such modes would have to start out at a larger amplitude in order for them 
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Run 


(-t>x; J^y, J^z) 


{N. 


^.Ny 


,iV.) 


^d/^g 


V^Vkq/cso 


Hd/Hg 


Initial Rimin 


-^rras 


Final Rimin 


01 


(0.1,0.4,0.4) 


(32,128,256) 


0.01 


0.1 


0.01 


0.0675 


0.01 


0.42 


02 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.1 


0.01 


0.0675 


0.004 


Stable 


03 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.1 


0.01 


0.0675 


0.001 


Stable 


04 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.1 


0.02 


0.270 


0.1 


Stable 


05 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.1 


0.02 


0.270 


0.2 


Stable 


06 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.1 


0.015 


0.152 


0.1 


0.23 


07 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.1 


0.015 


0.152 


0.04 


Stable 


08 


(0.05,0.2,0.2) 


(32 


128 


256) 


0.01 


0.1 


0.005 


0.0169 


10-6 


0.28 


09 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.1 


0.005 


0.0169 


10-6 


0.38 


10 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.02 


0.1 


0.01 


0.0675 


0.1 


0.36 


11 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.04 


0.1 


0.01 


0.0675 


0.1 


0.43 


12 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.08 


0.1 


0.01 


0.0675 


0.1 


0.47 


13 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.02 


0.1 


0.02 


0.270 


0.1 


Stable 


14 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.04 


0.1 


0.02 


0.270 


0.1 


Stable 


15 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.08 


0.1 


0.02 


0.270 


0.1 


Stable 


16 


(0.05,0.2,0.2) 


(32 


128 


256) 


0.02 


0.1 


0.005 


0.0169 


0.001 


0.28 


17 


(0.05,0.2,0.2) 


(32 


128 


256) 


0.04 


0.1 


0.005 


0.0169 


0.001 


0.30 


18 


(0.05,0.2,0.2) 


(32 


128 


256) 


0.08 


0.1 


0.005 


0.0169 


0.001 


0.40 


19 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.2 


0.02 


0.0675 


0.01 


0.21 


20 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.2 


0.04 


0.313 


0.1 


Stable 


21 


(0.1,0.4,0.4) 


(32 


128 


256) 


0.01 


0.2 


0.01 


0.0169 


0.01 


0.22 


22 


(0.05,0.2,0.2) 


(32 


128 


256) 


0.01 


0.05 


0.01 


0.270 


0.04 


Stable 


23 


(0.05,0.2,0.2) 


(32 


128 


256) 


0.01 


0.05 


0.005 


0.0675 


0.04 


0.32 


24 


(0.05,0.2,0.2) 


(32 


128 


256) 


0.01 


0.05 


0.0025 


0.0169 


0.04 


0.32 


25 


(0.1,0.4,0.4) 


(64 


256 


512) 


0.01 


0.1 


0.01 


0.0675 


0.01 


0.22 


26 


(0.1,0.4,0.4) 


(64 


256 


512) 


0.01 


0.1 


0.01 


0.0675 


0.001 


Stable 


27 


(0.05,0.2,0.2) 


(64 


256 


512) 


0.01 


0.1 


0.005 


0.0169 


10-6 


0.21 



Table 1: 3D Simulations of settled dust layers. 



-24- 



to grow to sufficient amplitude to trigger nonlinear effects before the shear damped further 
growth. This is indeed the case as the critical amplitude is roughly an order of magnitude 
higher than the for the layer half as thick. 

Runs 04 and 05 are for layers that are twice as thick as in Run 25: Hd/^g = 0.01, 
VoVkq/csq = 0.1, Hd/Hg = 0.02, /ij = 0.5, and Rimin = 0.270. The amplitudes of initial 
perturbations were Arms = 0.1 and Arms = 0.2, respectively. These layers are close to the 
critical Richardson number for stability in the absence of Coriolis force or shear, but would be 
unstable with the Coriolis force and no horizontal shear. In 3D simulations with horizontal 
shear, these layers are found to be stable to even relatively large amplitude perturbations. 
Thus, it appear s that the high- Richardson -number unstable flows with the Coriolis force ffist 



investigated by iGomez &: Ostrikerl (120051 ) are stabilized by the horizontal shear. 



For thinner layers than those in Run 25, however, the amplitude threshold practically 
vanishes. Runs 08, 09, and 27 are for a layer initially half as thick (Richardson number four 
times smaller): Hd/Hg = 0.01, r/oV/^o/cso = 0.1, Hd/Hg = 0.005, /ij = 2.0, and Rimin = 
0.0169 . The amplitude of perturbations was only 10~^, yet the layers were still unstable. 
Figure [H] shows the time evolution of this layer. There may indeed be a threshold, but it would 
be so low as to be practically irrelevant to the evolution of dust layers in real protoplanetary 
disk environments. 



4.2. Non-linear mixing and final states 

In simulations without horizontal shear, the turbulence ffiled the computational domain 
and almost completely mixed the dust with the gas. There was little evidence of a remaining 
dust layer and the dust-to-gas ratio was nearly uniform. However, in 3D simulations with 
horizontal shear, the turbulence and mixing were locally confined to a region above and 
below the original layer, resulting in the formation of a new, thicker layer. Figure [TOh shows 
the horizontally-averaged dust-to-gas ratio as a function of height for the initial and final 
dust layers in Runs 25 and 27. It is interesting to note that the final layers have very nearly 
the same profile even though they started out with different initial widths with very different 
growth rates. In Figure [TOb . the Richardson number as a function of height is graphed for 
these same two runs. Because the Richardson number involves the ratio of derivatives that 
both vanish at the midplane, it is not clear that the few points that have very low Richardson 
number near the midplane are not just numerical outliers. If we ignore those few points, then 
it appears that both of these runs result in minimum Richardson numbers right around the 
canonical value of one-quarter. The final column of Table [1] shows the numerically computed 
minimum Richardson number for the cases where the dust layer was unstable. The lower 
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Fig. 8. — 3D nonlinear evolution of Kelvin- Helmholtz instability with Coriolis force and 
horizontal shear. In this simulation, S^/S^ = 0.01 and fjoVKo/cso = 0.1. The initial dust 
scale height is Hd/Hg = 0.01, corresponding to a peak local dust-to-gas ratio in the midplane 
of Hq = 1.0 and a minimum Richardson number of Rimin = 0.0675. The amplitude of initial 
perturbations was Arms = 10^^. The first column illustrates the evolution of the local dust- 
to-gas ratio /i (deep red = 1.0, deep blue = 0); the second column shows the evolution of the 
radial component of vorticity u^ (red = vorticity that points into the page, blue = vorticity 
that points out of the page). The time interval between frames is 3.4 orbital periods. 
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Fig. 9. — 3D nonlinear evolution of Kelvin- Helmholtz instability with Coriolis force and 
horizontal shear. In this simulation, S^/S^ = 0.01 and tiqVkq/cso = 0.1. The initial dust 
scale height is Hd/Hg = 0.005, corresponding to a peak local dust-to-gas ratio in the midplane 
of Hq = 2.0 and a minimum Richardson number of Rimm = 0.0169. The amplitude of initial 
perturbations was Arms = 10^^. The first column illustrates the evolution of the local dust- 
to-gas ratio /i (deep red = 2.0, deep blue = 0); the second column shows the evolution of the 
radial component of vorticity u^ (red = vorticity that points into the page, blue = vorticity 
that points out of the page). The time interval between frames is 3.4 orbital periods. 
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resolution runs usually result in layers with minimum Richardson numbers that are slightly 
larger, in the 0.35 - 0.40 range. If the resulting turbulence yields a new layer that is thicker 
than the critical thickness, the dust has no way to re-sediment because these simulations 
are in the limit of perfect dust-to-gas coupling. In Figure [TT], the final horizontally-averaged 
dust-to-gas ratio for all the runs with unstable dust layers is plotted as a function of height. 
The axes are scaled in such a way so that profiles with the same minimum Richardson 
number coincide. The runs included in this figure include simulations with different global 
dust-to-gas ratios (Runs 10-18), different global gas radial pressure gradients (Runs 19-24), 
different layer widths, and different amounts of initial perturbations. Surprisingly, the vast 
majority of these simulations resulted in final dust layers with nearly the same minimum 
Richardson number. 



5. DISCUSSION & FUTURE WORK 

We revisited the case of no horizontal shear, both with and without the Coriolis force. 
The case without the Coriolis force is the most similar to classic Kelvin- Helmholtz instability, 
and the critical Richardson number for the onset of instability is close to the expected 
value of one-quarter, consistent with there being sufficient kinetic energy in the shear to 
lift t he heavier fluid out of the gravitational well and mix it w i th th e overlaying lighter 



fluid (JGaraud fc Linl 120041 ) . As first noted by iGomez fc Ostrikerl (120051 ) , the case with the 



Coriolis force is surprisingly different, with instability occurring at much higher Richardson 
numbers. However, the wavelengths of the most u nstable eigenmodes for t hese thicker layers 



is smaller than the thickness of the layer itself. Ilshitsu fc Sekiyal (120031 ) showed that the 
horizontal shear is able to eventually stabilize unstable eigenmodes by shearing them out to 
high wavenumber. However, unstable eigenmodes are able to grow for a period of time; the 
question is whether they can grow to a large enough amplitude to trigger nonlinear effects 
and disrupt the dust layer before they are damped as they are sheared. 

In order to investigate the competing influences of the Coriolis force (destabilizing) 
and horizontal shear (sta bilizing), we used a 3D spectral, anelastic, shearing-box code 



( jBarranco fc MarcusI 120061 ) to simulate settled dust layers in the limit of perfect dust-gas 
coupling. We find that the stability of dust layers depends on the amplitude of initial per- 
turbations: small perturbations grow for a period of time, but are damped before they reach 
sufficient magnitude to trigger nonlinear effects; whereas larger amplitude perturbations are 
able to grow to magnitudes that can disrupt the dust layer, resulting in turbulence and mix- 
ing. Kelvin-Helmholtz instability in thicker layers has slower growth rates, implying that 
the magnitude of initial amplitudes would have to be larger. This was seen in Runs 06 
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and 07 in which the critical threshold for perturbations was an order of magnitude larger 
for a layer that was 50% thicker. Thick dust layers that are stable according to the classic 
Richardson number criterion, but are unstable with the addition of the Coriolis force exhibit 
Kelvin-Helmholtz instability with very slow growth rates and at high spatial wavenumber 
(see Figure [5]). Three-dimensional simulations of these layers indicate that the horizontal 
shear is able to damp the instability no matter how large the initial perturbations for these 
layers (see Runs 04 and 05). Thinner layers, on the other hand, have very fast growth rates, 
and the threshold amplitude is so low as to be practically irrelevant in real protoplanetary 
disks where there is no doubt fluctuations of such low magnitude. 

In two-dimensional simulations of unstable dust layers, the turbulence that develops 
fills the computational domain and leads to large-scale mixing. Characteristic of "inverse 
cascades" in two-dimensional turbulence, small eddies merge with other eddies to form larger 
coherent vortices which chaotically advect the dust, homogenizing the dust-to-gas ratio. In 
three-dimensional simulations, the turbulence is locally confined to a region right around 
the original unstable dust layer and does not propagate through the rest of the computa- 
tional domain. The unstable thin layers evolve to thicker, more stable layers with minimum 
Richardson numbers tantalizingly close to the value of one-quarter for the onset of instability 
in the absence of the Coriolis force and horizontal shear. These results hold when the global 
dust to gas ratio S^/S^ and the global gas radial pressure gradient r/oV^-o/cso are varied. 
There is no reason to expect that the unstable layers in these simulations would evolve to a 
final state that just happens to be at the critical thickness for stability. One could imagine 
that the turbulence is so efficient at mixing that it results in layers whose widths are signif- 
icantly thicker than thinnest stable layer. Because the simulations presented here are in the 
limit of perfect dust-to-gas coupling, no further sedimentation of the dust is allowed so that 
the thicker final states cannot settle to the critical state. 

Future work will involve two significant improvements. First, we will relax the perfect 
dust-to-gas coupling assumption and allow there to be a finite v alue for the stoppi ng time. 
Dust will be treated as a second fiuid with its own velocity field (ICuzzi et al.lll993l ). Layers 
that are unstable will devel op turbulence and re-m ix the dust with the gas, but then would 
be allowed to further settle. iJohansen et al.l (120061 ) performed 2D simulations with two fiuids 
and showed that layers evolved toward a self-regulated state in which further settling was 
inhibited by turbulence generated by Kelvin-Helmholtz instability, maintaining the layer in 
a dynamic equilibrium right at the critical thickness for instability. However, because their 
simulations were 2D, they found the dust layers had very high Richardson number similar to 
the layers investigated by iGomez fc Ostriken (120051 ) . Second, we wi ll add the effects of non- 
ideal magnetohydrodynamics with finite resistance. [Turner et al.l (120071 ) investigated the 
magneto-rotational instability in protoplanetary disks with "dead zones" in the midplane 
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where the ionization is too low to couple the fluid to the magnetic fields. However, they find 
that turbulence originating in the cosmic-ray-ionized surface layers can mix free charges into 
the interior and weakly couple the midplane gas to the magnetic fields, effectively eliminating 
the dead zone. Their analysis does not include the role that dust grains play in removing 
free charges and reducing the ionization. If turbulence lofts the particles throughout the 
disk, ionization can be suppressed, which will have the tendency to decouple the interior 
gas from the magnetic fields, allowing the dead zone to reform. Dust particle could then 
re-settle, ionization could increase, re-coupling the gas to the fields and generating a new 
phase of turbulence. We plan to investigate if there is indeed a limit cycle to the formation 
and destruction of the dead zone in protoplanetary disks. 
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Fig. 10. — Horizontally-averaged dust-to-gas ratio fi and Richardson number Ri as a function 
of height z for initial and final dust layers for runs 25 and 27. 
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Fig. 11. — Horizontally-averaged dust-to-gas ratio as a function of height for runs with 
unstable dust layers in Table [H Axes are scaled in such a way so that layers that have the 
same Richardson number would lay on top of one another. 



